#include "NoisePerlin.hpp"

namespace zzz{
template<>
const Vector<3,float> NoisePerlin<float> ::vector_[16] =
{
    Vector<3,float>(1, 1, 0),
    Vector<3,float>(-1, 1, 0),
    Vector<3,float>(1,-1, 0),
    Vector<3,float>(-1,-1, 0),
    Vector<3,float>(1, 0, 1),
    Vector<3,float>(-1, 0, 1),
    Vector<3,float>(1, 0,-1),
    Vector<3,float>(-1, 0,-1),
    Vector<3,float>(0, 1, 1),
    Vector<3,float>(0,-1, 1),
    Vector<3,float>(0, 1,-1),
    Vector<3,float>(0,-1,-1),
    Vector<3,float>(1, 1, 0),
    Vector<3,float>(-1, 1, 0),
    Vector<3,float>(0,-1, 1),
    Vector<3,float>(0,-1,-1)
};

template<>
float NoisePerlin<float>::PerlinNoise(const zuint seed,const zuint x,const zuint fractionalBits,const float fx)
{
  // Some sanity checking in debug mode
  assert(fx>=0.0 && fx<=1.0);
  assert(fractionalBits<32);
  // Determine the integer co-ordinates of the lattice
  // interval to be interpolated.  Note that the lattice
  // tiles seamlessly and fractionalBits are discarded.
  const zuint px[2] = { x>>fractionalBits, (x+(1<<fractionalBits))>>fractionalBits };
  // The floating-point fx displacement is
  // scaled down to be less significant than the
  // fractional components of x,y and z
  const float subFraction = 1/float(1<<fractionalBits);
  // Combine the fractional bits of x and y
  // with the floating-point displacements
  // to determine the location within the square
  float fx0;
  Fraction(fx0,(x<<(32-fractionalBits)));
  fx0 += fx*subFraction;
  // Determine the complement value of fx0
  const float fx1 = fx0-1;
  // Apply ease-in ease-out curve, rather
  // than interpolating linearly
  const float wx = Scurve(fx0);
  // Temporary storage for interpolation
  zuint i;
  float vx0, vx1;
  // Interpolate along edge from origin in +x direction
  i = Hash(seed,px[0]);
  vx0 = vector_[i].x()*fx0;
  i = Hash(seed,px[1]);
  vx1 = vector_[i].x()*fx1;
  return vx0 + wx*(vx1-vx0);
}

template<>
float NoisePerlin<float>::PerlinNoise(const zuint seed,const zuint x,const zuint y,const zuint fractionalBits,const float fx,const float fy)
{
  // Some sanity checking in debug mode
  assert(fx>=0.0 && fx<=1.0);
  assert(fy>=0.0 && fy<=1.0);
  assert(fractionalBits<32);
  // Determine the integer co-ordinates of the lattice
  // square to be interpolated.  Note that the lattice
  // tiles seamlessly and fractionalBits are discarded.
  const zuint px[2] = { x>>fractionalBits, (x+(1<<fractionalBits))>>fractionalBits };
  const zuint py[2] = { y>>fractionalBits, (y+(1<<fractionalBits))>>fractionalBits };
  // The floating-point fx and fy displacements are
  // scaled down to be less significant than the
  // fractional components of x,y and z
  const float subFraction = 1/float(1<<fractionalBits);
  // Combine the fractional bits of x and y
  // with the floating-point displacements
  // to determine the location within the square
  float fx0;
  Fraction(fx0,(x<<(32-fractionalBits)));
  fx0 += fx*subFraction;
  float fy0;
  Fraction(fy0,(y<<(32-fractionalBits)));
  fy0 += fy*subFraction;
  // Determine the complement values of fx0 and fy0
  const float fx1 = fx0-1;
  const float fy1 = fy0-1;
  // Apply ease-in ease-out curve, rather
  // than interpolating linearly
  const float wx = Scurve(fx0);
  const float wy = Scurve(fy0);
  // Temporary storage for interpolation
  zuint i;
  float vx0, vx1, vy0, vy1;
  // Interpolate along edge from origin in +x direction
  i = Hash(seed,px[0],py[0]);
  vx0 = vector_[i].x()*fx0 + vector_[i].y()*fy0;
  i = Hash(seed,px[1],py[0]);
  vx1 = vector_[i].x()*fx1 + vector_[i].y()*fy0;
  vy0 = vx0 + wx*(vx1-vx0);
  // Interpolate along edge from +y in +x direction
  i = Hash(seed,px[0],py[1]);
  vx0 = vector_[i].x()*fx0 + vector_[i].y()*fy1;
  i = Hash(seed,px[1],py[1]);
  vx1 = vector_[i].x()*fx1 + vector_[i].y()*fy1;
  vy1 = vx0 + wx*(vx1-vx0);
  // Interpolate between the two edges
  return vy0 + wy*(vy1-vy0);
}

template<>
float NoisePerlin<float>::PerlinNoise(const zuint seed,const zuint x,const zuint y,const zuint z,const zuint fractionalBits,const float fx,const float fy,const float fz)
{
  // Some sanity checking in debug mode
  assert(fx>=0.0 && fx<=1.0);
  assert(fy>=0.0 && fy<=1.0);
  assert(fz>=0.0 && fz<=1.0);
  assert(fractionalBits<32);
  // Determine the integer co-ordinates of the lattice
  // cube to be interpolated.  Note that the lattice
  // tiles seamlessly and fractionalBits are discarded.
  const zuint px[2] = { x>>fractionalBits, (x+(1<<fractionalBits))>>fractionalBits };
  const zuint py[2] = { y>>fractionalBits, (y+(1<<fractionalBits))>>fractionalBits };
  const zuint pz[2] = { z>>fractionalBits, (z+(1<<fractionalBits))>>fractionalBits };
  // The floating-point fx,fy,fz displacements are
  // scaled down to be less significant than the
  // fractional components of x,y and z
  const float subFraction = 1/float(1<<fractionalBits);
  // Combine the fractional bits of x,y and z
  // with the floating-point displacements
  // to determine the location within the cube
  float fx0;
  Fraction(fx0,(x<<(32-fractionalBits)));
  fx0 += fx*subFraction;
  float fy0;
  Fraction(fy0,(y<<(32-fractionalBits)));
  fy0 += fy*subFraction;
  float fz0;
  Fraction(fz0,(z<<(32-fractionalBits)));
  fz0 += fz*subFraction;
  // Determine the complement values of fx0,fy0 and fz0
  const float fx1 = fx0-1;
  const float fy1 = fy0-1;
  const float fz1 = fz0-1;
  // Apply ease-in ease-out curve, rather
  // than interpolating linearly
  const float wx = Scurve(fx0);
  const float wy = Scurve(fy0);
  const float wz = Scurve(fz0);
  // Temporary storage for interpolation
  zuint i;
  float vx0, vx1, vy0, vy1, vz0, vz1;
  // Interpolate along edge from origin in +x direction
  i = Hash(seed,px[0],py[0],pz[0]);
  vx0 = vector_[i].x()*fx0 + vector_[i].y()*fy0 + vector_[i].z()*fz0;
  i = Hash(seed,px[1],py[0],pz[0]);
  vx1 = vector_[i].x()*fx1 + vector_[i].y()*fy0 + vector_[i].z()*fz0;
  vy0 = vx0 + wx*(vx1-vx0);
  // Interpolate along edge from +y in +x direction
  i = Hash(seed,px[0],py[1],pz[0]);
  vx0 = vector_[i].x()*fx0 + vector_[i].y()*fy1 + vector_[i].z()*fz0;
  i = Hash(seed,px[1],py[1],pz[0]);
  vx1 = vector_[i].x()*fx1 + vector_[i].y()*fy1 + vector_[i].z()*fz0;
  vy1 = vx0 + wx*(vx1-vx0);
  // Interpolate between the two edges
  vz0 = vy0 + wy*(vy1-vy0);
  // Interpolate along edge from +z in +x direction
  i = Hash(seed,px[0],py[0],pz[1]);
  vx0 = vector_[i].x()*fx0 + vector_[i].y()*fy0 + vector_[i].z()*fz1;
  i = Hash(seed,px[1],py[0],pz[1]);
  vx1 = vector_[i].x()*fx1 + vector_[i].y()*fy0 + vector_[i].z()*fz1;
  vy0 = vx0 + wx*(vx1-vx0);
  // Interpolate along edge from +y+z in +x direction
  i = Hash(seed,px[0],py[1],pz[1]);
  vx0 = vector_[i].x()*fx0 + vector_[i].y()*fy1 + vector_[i].z()*fz1;
  i = Hash(seed,px[1],py[1],pz[1]);
  vx1 = vector_[i].x()*fx1 + vector_[i].y()*fy1 + vector_[i].z()*fz1;
  vy1 = vx0 + wx*(vx1-vx0);
  // Interpolate between the two +z edges
  vz1 = vy0 + wy*(vy1-vy0);
  // Finally, interpolate betwen the two planes
  return vz0 + wz*(vz1-vz0);
}
//////////////////////////////////////////////////////////////////////////
template<>
const Vector<3,double> NoisePerlin<double>::vector_[16] =
{
  Vector<3,double>(1, 1, 0),
  Vector<3,double>(-1, 1, 0),
  Vector<3,double>(1,-1, 0),
  Vector<3,double>(-1,-1, 0),
  Vector<3,double>(1, 0, 1),
  Vector<3,double>(-1, 0, 1),
  Vector<3,double>(1, 0,-1),
  Vector<3,double>(-1, 0,-1),
  Vector<3,double>(0, 1, 1),
  Vector<3,double>(0,-1, 1),
  Vector<3,double>(0, 1,-1),
  Vector<3,double>(0,-1,-1),
  Vector<3,double>(1, 1, 0),
  Vector<3,double>(-1, 1, 0),
  Vector<3,double>(0,-1, 1),
  Vector<3,double>(0,-1,-1)
};

template<>
double NoisePerlin<double>::PerlinNoise(const zuint seed,const zuint x,const zuint fractionalBits,const double fx)
{
  // Some sanity checking in debug mode
  assert(fx>=0.0 && fx<=1.0);
  assert(fractionalBits<32);
  // Determine the integer co-ordinates of the lattice
  // interval to be interpolated.  Note that the lattice
  // tiles seamlessly and fractionalBits are discarded.
  const zuint px[2] = { x>>fractionalBits, (x+(1<<fractionalBits))>>fractionalBits };
  // The floating-point fx displacement is
  // scaled down to be less significant than the
  // fractional components of x,y and z
  const double subFraction = 1/double(1<<fractionalBits);
  // Combine the fractional bits of x and y
  // with the floating-point displacements
  // to determine the location within the square
  double fx0;
  Fraction(fx0,(x<<(32-fractionalBits)));
  fx0 += fx*subFraction;
  // Determine the complement value of fx0
  const double fx1 = fx0-1;
  // Apply ease-in ease-out curve, rather
  // than interpolating linearly
  const double wx = Scurve(fx0);
  // Temporary storage for interpolation
  zuint i;
  double vx0, vx1;
  // Interpolate along edge from origin in +x direction
  i = Hash(seed,px[0]);
  vx0 = vector_[i].x()*fx0;
  i = Hash(seed,px[1]);
  vx1 = vector_[i].x()*fx1;
  return vx0 + wx*(vx1-vx0);
}

template<>
double NoisePerlin<double>::PerlinNoise(const zuint seed,const zuint x,const zuint y,const zuint fractionalBits,const double fx,const double fy)
{
  // Some sanity checking in debug mode
  assert(fx>=0.0 && fx<=1.0);
  assert(fy>=0.0 && fy<=1.0);
  assert(fractionalBits<32);
  // Determine the integer co-ordinates of the lattice
  // square to be interpolated.  Note that the lattice
  // tiles seamlessly and fractionalBits are discarded.
  const zuint px[2] = { x>>fractionalBits, (x+(1<<fractionalBits))>>fractionalBits };
  const zuint py[2] = { y>>fractionalBits, (y+(1<<fractionalBits))>>fractionalBits };
  // The floating-point fx and fy displacements are
  // scaled down to be less significant than the
  // fractional components of x,y and z
  const double subFraction = 1/double(1<<fractionalBits);
  // Combine the fractional bits of x and y
  // with the floating-point displacements
  // to determine the location within the square
  double fx0;
  Fraction(fx0,(x<<(32-fractionalBits)));
  fx0 += fx*subFraction;
  double fy0;
  Fraction(fy0,(y<<(32-fractionalBits)));
  fy0 += fy*subFraction;
  // Determine the complement values of fx0 and fy0
  const double fx1 = fx0-1;
  const double fy1 = fy0-1;
  // Apply ease-in ease-out curve, rather
  // than interpolating linearly
  const double wx = Scurve(fx0);
  const double wy = Scurve(fy0);
  // Temporary storage for interpolation
  zuint i;
  double vx0, vx1, vy0, vy1;
  // Interpolate along edge from origin in +x direction
  i = Hash(seed,px[0],py[0]);
  vx0 = vector_[i].x()*fx0 + vector_[i].y()*fy0;
  i = Hash(seed,px[1],py[0]);
  vx1 = vector_[i].x()*fx1 + vector_[i].y()*fy0;
  vy0 = vx0 + wx*(vx1-vx0);
  // Interpolate along edge from +y in +x direction
  i = Hash(seed,px[0],py[1]);
  vx0 = vector_[i].x()*fx0 + vector_[i].y()*fy1;
  i = Hash(seed,px[1],py[1]);
  vx1 = vector_[i].x()*fx1 + vector_[i].y()*fy1;
  vy1 = vx0 + wx*(vx1-vx0);
  // Interpolate between the two edges
  return vy0 + wy*(vy1-vy0);
}

template<>
double NoisePerlin<double>::PerlinNoise(const zuint seed,const zuint x,const zuint y,const zuint z,const zuint fractionalBits,const double fx,const double fy,const double fz)
{
  // Some sanity checking in debug mode
  assert(fx>=0.0 && fx<=1.0);
  assert(fy>=0.0 && fy<=1.0);
  assert(fz>=0.0 && fz<=1.0);
  assert(fractionalBits<32);
  // Determine the integer co-ordinates of the lattice
  // cube to be interpolated.  Note that the lattice
  // tiles seamlessly and fractionalBits are discarded.
  const zuint px[2] = { x>>fractionalBits, (x+(1<<fractionalBits))>>fractionalBits };
  const zuint py[2] = { y>>fractionalBits, (y+(1<<fractionalBits))>>fractionalBits };
  const zuint pz[2] = { z>>fractionalBits, (z+(1<<fractionalBits))>>fractionalBits };
  // The floating-point fx,fy,fz displacements are
  // scaled down to be less significant than the
  // fractional components of x,y and z
  const double subFraction = 1/double(1<<fractionalBits);
  // Combine the fractional bits of x,y and z
  // with the floating-point displacements
  // to determine the location within the cube
  double fx0;
  Fraction(fx0,(x<<(32-fractionalBits)));
  fx0 += fx*subFraction;
  double fy0;
  Fraction(fy0,(y<<(32-fractionalBits)));
  fy0 += fy*subFraction;
  double fz0;
  Fraction(fz0,(z<<(32-fractionalBits)));
  fz0 += fz*subFraction;
  // Determine the complement values of fx0,fy0 and fz0
  const double fx1 = fx0-1;
  const double fy1 = fy0-1;
  const double fz1 = fz0-1;
  // Apply ease-in ease-out curve, rather
  // than interpolating linearly
  const double wx = Scurve(fx0);
  const double wy = Scurve(fy0);
  const double wz = Scurve(fz0);
  // Temporary storage for interpolation
  zuint i;
  double vx0, vx1, vy0, vy1, vz0, vz1;
  // Interpolate along edge from origin in +x direction
  i = Hash(seed,px[0],py[0],pz[0]);
  vx0 = vector_[i].x()*fx0 + vector_[i].y()*fy0 + vector_[i].z()*fz0;
  i = Hash(seed,px[1],py[0],pz[0]);
  vx1 = vector_[i].x()*fx1 + vector_[i].y()*fy0 + vector_[i].z()*fz0;
  vy0 = vx0 + wx*(vx1-vx0);
  // Interpolate along edge from +y in +x direction
  i = Hash(seed,px[0],py[1],pz[0]);
  vx0 = vector_[i].x()*fx0 + vector_[i].y()*fy1 + vector_[i].z()*fz0;
  i = Hash(seed,px[1],py[1],pz[0]);
  vx1 = vector_[i].x()*fx1 + vector_[i].y()*fy1 + vector_[i].z()*fz0;
  vy1 = vx0 + wx*(vx1-vx0);
  // Interpolate between the two edges
  vz0 = vy0 + wy*(vy1-vy0);
  // Interpolate along edge from +z in +x direction
  i = Hash(seed,px[0],py[0],pz[1]);
  vx0 = vector_[i].x()*fx0 + vector_[i].y()*fy0 + vector_[i].z()*fz1;
  i = Hash(seed,px[1],py[0],pz[1]);
  vx1 = vector_[i].x()*fx1 + vector_[i].y()*fy0 + vector_[i].z()*fz1;
  vy0 = vx0 + wx*(vx1-vx0);
  // Interpolate along edge from +y+z in +x direction
  i = Hash(seed,px[0],py[1],pz[1]);
  vx0 = vector_[i].x()*fx0 + vector_[i].y()*fy1 + vector_[i].z()*fz1;
  i = Hash(seed,px[1],py[1],pz[1]);
  vx1 = vector_[i].x()*fx1 + vector_[i].y()*fy1 + vector_[i].z()*fz1;
  vy1 = vx0 + wx*(vx1-vx0);
  // Interpolate between the two +z edges
  vz1 = vy0 + wy*(vy1-vy0);
  // Finally, interpolate betwen the two planes
  return vz0 + wz*(vz1-vz0);
}
}
